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Abstract 

Microbial populations in the natural environment are likely to experience growth 
conditions very different from those of a typical laboratory experiment. In par- 
ticular, removal rates of biomass and substrate are unlikely to be balanced 
under realistic environmental conditions. Here, we consider a single population 
growing on a substrate under conditions where the removal rates of substrate 
and biomass are not necessarily equal. For a large population, with determin- 
istic growth dynamics, our model predicts that this system can show transient 
(damped) oscillations. For a small population, demographic noise causes these 
oscillations to be sustained indefinitely. These oscillations arise when the dy- 
namics of changes in biomass are faster than the dynamics of the substrate, for 
example, due to a high microbial death rate and/or low substrate flow rates. 
We show that the same mechanism can produce sustained stochastic oscillations 
in a two-species, nutrient-cycling microbial ecosystem. Our results suggest that 
oscillatory population dynamics may be a common feature of small microbial 
populations in the natural environment, even in the absence of complex inter- 
species interactions. 

Keywords: microbial growth, chemostat, demographic fluctuations, noise, 
microbial ecology, power spectrum 



1. Introduction 

Microbial populations play an essential role in many important processes 
in the natural and human environments, including nitrification of soil, carbon 
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processing in ocean food chains and wastewater treatment. Extensive stud- 
ies have informed our understanding of microbial population dynamics under 
standard conditions in the laboratory, yet outside the lab, different conditions 
may hold. In particular, nutrient supply may be unpredictable, nutrient and 
biomass removal rates may not balance, and population sizes may be small 
so that random fluctuations due to birth and death events play an important 
role. Recent studies have shown that microbial communities can undergo un- 
predictable divergence from similar initial conditions [l|, dramatic fluctuations 
in species composition f5| and even chaotic dynamics [^Bl ■ Stochastic models 
for microbial ecosystem dynamics will be needed to obtain detailed understand- 
ing of such results. Here, we take a first step towards developing such models, 
by studying the effects of unbalanced biomass and substrate fiow rates, and of 
demographic noise, on microbial population dynamics. 

Microbial population dynamics are often studied in the context of the chemo- 
stat 0, @]: a well-stirred vessel in which a population is maintained under 
steady-state conditions, with constant inflow of substrate at fixed concentration, 
balanced by constant outfiow of the vessel's contents (biomass and substrate). 
Assuming the standard Monod relation between growth rate and substrate con- 
centration, the deterministic equations for biomass and substrate concentrations 
in the chemostat do not show either transient or sustained oscillations [8]. A 
large body of theoretical work has shown that oscillations can occur in extended 
versions of the classic chemostat equations, which include factors such as time 
delays [ 9 | . [l0l , feedback control [III , growth inhibition by substrate or prod- 
ucts [13I Il4|. age structure within the population fis'l, periodic nutrient input 
or washout rates ^M^, and predator/prey interactions among populations 



within the chemostat |19| - and indeed, oscillations have been observed in a 
number of chemostat experiments However, almost all these models 

assume equal rates of substrate and biomass removal from the system (which is 
a consequence of the standard chemostat setup) , and deterministic population 
dynamics. In this paper, we take a different approach: we use a very simple 
growth model, but take into account both imbalances in substrate and biomass 
flow rates, and fluctuations due to small population size. 

For small microbial populations, individual birth and death events give rise 
to demographic fluctuations which may be significant compared to the total 
population size. For natural communities growing in microenvironments such 
as the interstices between soil grains or the surfaces of ocean particles, these de- 
mographic fluctuations may well be relevant. Moreover, recent work has shown 
that, in populations which expand to colonize new spatial territories, demo- 
graphic fluctuations, due to small numbers at the boundary, can have dramatic 
consequences for the genetic structure of the population, even for large popula- 
tions H, 23] • From a theoretical point of view, it is well known that fluctuations 
can have important qualitative effects on the behaviour of dynamical systems 
25- 23|- In particular, Newman and McKane have shown recently that sustained 



oscillations can arise in stochastic dynamical systems driven by intrinsic noise, 
whose corresponding deterministic dynamical equations lead only to transient 
(damped) oscillations [27-29]. This effect has been observed in a range of model 
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systems including predator-prey models 27h29||. models for infection dynamics 
|31| , cooperative gam es l 32l . |33| and chemical reaction models for genetic and 



metabolic regulation l34l [35| . These oscillations arise because intrinsic noise ex- 
cites the underlying damped oscillatory degrees of freedom of the system; they 
are self-sustaining, with all initial conditions leading to the same ensemble of 
stochastic orbits, and (similar to a limit-cycle) a characteristic amplitude and 
frequency are set by the parameters of the system. In this paper, we show that 
the same mechanism can lead to sustained stochastic oscillations for a simple, 
one-population microbial growth model with unbalanced rates of substrate and 
biomass removal, and we suggest that it is likely to be a generic feature of 
microbial ecosystems with small population sizes □. 

In Section [5J we analyse a deterministic model for a single microbial popu- 
lation growing on a substrate, with unbalanced rates of biomass and substrate 
removal; we show that this system can show transient oscillations. In Section 
[3l we show that demographic noise can cause sustained stochastic oscillations 
in this system. In Section |4l we extend our analysis to a simple two-species, 
nutrient-cycling model. Finally, we present our conclusions in Section O 



2. A single population: deterministic model 

We first consider, using deterministic equations, the dynamics of a single mi- 
crobial population with biomass concentration x{t) (units of microbial cells per 
litre), which consumes a substrate of concentration s{t) (units of /iM). Substrate 
flows into the system at a constant rate b (units of concentration/time) and is 
removed from the system with rate constant R (units of time^^), while biomass 
is removed with rate constant D (also with units of time^^). In a chemostat, the 
rates of substrate and biomass removal are equal {R = Z?), since the well- mixed 
contents of the chemostat are pumped out at a constant rate, and the microbial 
death rate is assumed to be negligible. In the natural environment, however, the 
situation is more complex: microbes may be subject to significant (e.g. phage- 
mediated) killing, substrate may be removed via consumption by competing 
organisms, and microbes may avoid being washed away by adhering to a sur- 
face. We do not therefore constrain the removal rates of substrate and biomass 
to be equal. Since in this section we use a deterministic approach, neglecting 
demographic fluctuations, our analysis is appropriate for a large population; in 
Section [3] we present contrasting results for a small population. 

The equations governing the dynamics of our system are 

i — f{x,s)—^{s)x — Dx (1) 
s ~ g{x, s) — —j^{s)x + b — Rs (2) 



dx 
'dt 
ds 

dt 



^We note that this mechanism, in which stochastic oscillations are generated by intrinsic 
noise in the system dynamics, is distinct from the phenomenon of coherence resonance [5^, 
in which oscillations are caused by noise in an external driving force. 
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where we assume that microbial growth is related to substrate concentration by 
the Monod function [SGtiSTj: 



V d J ^ ^^'^^ " kTS 

\ / growth J\ -t t, 

in which v is the maximal growth rate and K is the substrate concentration 
at which the growth rate is half-maximal. Thus the growth rate is linearly 
proportional to substrate concentration when s <C if, but saturates at high 
substrate concentration s ^ K. The parameter 7 is the number of substrate 
molecules that need to be consumed to make one microbe (the inverse of the 
yield coefhcient (38|). 

Eqs (dl and ^ have a single non-trivial fixed point (for which i = s = 0) 

at 

(4) 



(v/D - 1) 
(b-Rs*) 



Starting from an arbitrary (nonzero) initial condition, the system will evolve 
towards this fixed point. To determine whether this happens monotonically or 
in an oscillatory manner, we analyse the system's dynamics close to the fixed 
point, by making the linear approximation (39| : 

6x \ _ f dxf dsf \ f Sx \ _ y ^ Sx \ 



Ss J \ dx9 dsg /(^.^.)t V '^^ / V 

where Sx = x — x* , 6s — s ~ s* , dzf and dsf are shorthand for df/dz and 
df/ds respectively, and J* — J(a;*,s*) denotes the Jacobian matrix of first- 
order partial derivatives, evaluated at the fixed point (a;*,s*). If Ui and U2 
are the eigenvectors of the Jacobian matrix J* at the fixed point, then Eq.(l5]) 
implies that the system evolves as ttie'^^* + 1126'^^*, where Ai and A2 are the 
corresponding eigenvalues of J*. Real and negative eigenvalues of J* indicate 
exponential relaxation to the fixed point, while complex eigenvalues with a 
negative real part indicate exponentially decaying oscillations as the system 
approaches its fixed point fsd]. Using Eqs ([T]), ([2|) and ([4]), we obtain: 



^( 1 JJL. ) (6) 



-7 -iP + x) 
where we have defined two dimensionless parameters: 

measures the rate of substrate removal relative to the rate of biomass removal, 
and 

-fx* ( d^x{s)\ 



D \ ds 
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measures the responsiveness of the microbial growth rate to changes in the 
substrate concentration. Increasing the substrate inflow rate b results in an 
increase in /3 (since x* depends on b). The eigenvalues A of this Jacobian are 
given by: 

A _ -(/3 + x)± V(/3 + x)^-4/? 
D~ 2 

The case x > 1 corresponds to a scenario where substrate is removed from 
the system faster than, or at the same rate as, biomass is removed. Under 
these circumstances, the eigenvalues A are real and negative for any value of 
/3, indicating that the system relaxes exponentially to its fixed point, with no 
oscillations (note that x — ^ corresponds to the chemostat case[§). However, 
if substrate removal is slower than biomass removal (x < 1), then the eigen- 
values A can be complex (with negative real part), implying that the system 
can undergo transient oscillations in biomass and substrate concentration as it 
approaches the fixed point. These oscillations occur over a range of parameter 
values corresponding to (/3 + x)^ < 4/3: i.e. for (3 in the range (5^ < (3 < /?+, 
where f3± = 2 — x ± 2-y/l — x- The coloured region in Figure [1] shows the region 
of /3 — X parameter space where transient oscillations are expected. 

To show that these transient oscillations are relevant for microbial popula- 
tions, we plot in Figure [2] simulated dynamical trajectories for biomass and sub- 
strate concentration, for a parameter set whose values are chosen to correspond 
approximately to Escherichia coli growing on glucose: v = lhr~^, D — 0.5hr~^, 
K = 1/zM and 7 = 1.8 x 10^" substrate molecules consumed to produce one 
microbe [s^l [S] [Si • Keeping these parameters fixed, the dimensionless param- 
eters /3 and x in Eqs. (O and (|8]) are controlled by the substrate flow rates 
with the following numerical values: (3 — b — R and x = 2i?; we choose {(3, x) 
combinations corresponding to the green circles in the phase diagram of Figure 
[TJ Significant transient oscillations are indeed observed in our simulations (on 
a timescale of typically tens of hours) for those parameter combinations which 
lie inside the coloured region of Figure [TJ 

Our analysis of the eigenvalues of the Jacobian (Eqn.®) also provides in- 
formation on the frequency of the transient oscillations that occur during 
relaxation to the fixed point. This is given by the imaginary part of the eigen- 
values A: 

^ - f v/4/3-(/3 + x)^- (10) 

The predicted frequency f2 is indicated by the contours in the phase diagram 
of Figure [TJ Increasing x - i-e. either increasing the rate of substrate removal 
or decreasing the rate of biomass removal - decreases the oscillation frequency, 
while the dependence on (3 is nonmonotonic, with a peak in VL for intermediate 
values of j3 (for fixed x)- 



•'if X = !> Eqns. [T] and [2] share a symmetry which reduces the number of dynamical degrees 
of freedom to one, precluding the possibility of oscillating solutions. 
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Figure 1: Summary of the predicted behaviour of the single population model as a function of 
the parameter combinations /3 and x- The coloured area shows the region of /3 — x parameter 
space where our deterministic model predicts transient oscillations during relaxation to the 
steady state (i.e. the region where the eigenvalues of the Jacobian are complex with negative 
real part). Within this region, the contours show the frequency of the damped oscillations 
(Eqn.l llOll '). The solid red line shows the boundary of the region where our stochastic model 
predicts sustained oscillations (i.e. the values of /3 and x for which there is a peak in the 
biomass power spectrum Eqn. JJTJ). The dotted line represents the contour for which the 
quality factor Q = 1 in Eqn.JTTJ. The green circles correspond to the parameter values used 
in the simulations of Figures[2l[l]and[6l where the logarithmic x-axis is schematically extended 
to zero in order to represent simulations for x = 0. 



To understand better the nature of the oscillations, we also compute the 
quality factor Q, which is a dimensionless measure of the extent to which oscil- 
lations are damped, and corresponds roughly to the number of oscillations that 
occur before the oscillations die away in the deterministic time series (Fig. [2]). 
For a simple harmonic oscillator, Q is the ratio of the energy stored to energy 
dissipated over an oscillation cycle. As we show in [Appendix A[ close to the 
fixed point, our deterministic model can be mapped onto the equations for a 
simple harmonic oscillator; this allows us to approximate the Q-factor of the 
oscillations as 

Q-^. (11) 

/3 + X 

Figure [3] shows Q (from Eq.lfTT])) as a function of /3, for several values of 
X- If /? > x^ the quality factor of the oscillations increases as /? decreases: 
reduced inflow of substrate (keeping other parameters fixed) will lead to more 
pronounced oscillations. However, if /3 < x, the opposite scenario holds; the 
quality factor of the oscillations will increase on increasing /? - i.e. on increasing 
the inflow of substrate, keeping all other parameters fixed. The dashed black line 
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Figure 2: Deterministic time series for the growth of Escherichia coli on glucose, calculated 
by numerical simulation of the coupled ODEs represented by Eqns.Q&lO, for parameters 
V = lhr~^, D = 0.5hr~^, K = 1/xM, 7 = 1.8 X 10^" substrate molecules. The 6 combinations 
of /3 and x correspond to the green circles in the /3 — x phase diagram of Figure [T] i.e. 
P = {0.01, 0.1, 1, 10} with X = (6 = {0.01, 0.1, 1, 10}/iMhr-l) and /3 = 0.1 with x = {0.25, 1} 
{b = 0.225/.tMhr-\ R = 0.125hr-l and b = OSfiMhr"^ , R = O.Shr'^). Note that in the 
legend, X = unless otherwise stated. 




Figure 3: Quality factor Q, predicted by Eo.l llll l. plotted as a function of /3 for various values 
of X- 



in the phase diagram of Figure [T] shows the contour of /3 and x for which Q = 1; 
inside this Une we expect to see that oscihations persist for a long time, while 
outside this line, oscillations, even though present, are strongly damped and 
rapidly decay. To understand this further, we rescale the substrate deviation 
from the fixed point as 6s' — ^^Ss and scale time by D, the death/removal 

rate of biomass. In these scaled units, the substrate relaxes like Ss' = —(/? + 
x){Sx + ds') and the biomass varies as Sx — -g^^Ss' . We can then identify the 
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relaxation rate constant of the substrate as = /? + x and that for the biomass 
as Xx = ^g+^j from which we see that 

A./A, = Q\ (12) 

In other words the more quickly the biomass responds to changes in substrate 
relative to the rate of relaxation of the substrate the higher the quality factor 
of the oscillations and the longer they will persist. Conversely, if the substrate 
relaxes more quickly than the biomass can respond, we see that the quality 
factor Q <C 1 and the oscillations will be very strongly damped (for those values 
of f3 and x expected to give rise to oscillations - as shown in Fig[T]). 

The underlying cause of the transient oscillations can be understood in terms 
of changes in balance between microbial growth and death and between con- 
sumption and net inflow of substrate. The phase plane plot of Fig. |3J^ shows 
how the biomass and substrate concentrations change during the approach to 
the fixed point, for trajectories starting from three different initial conditions, 
with j3 — 0.1 and x = 0. Transient oscillations are apparent from the fact that 
the trajectories spiral into the fixed point. The nullclines i — and s = 
(s = s* and b ~ Rs — -ffi{s)x respectively, obtained by setting Eqns.(IT]) and 
© to zero) , shown as dashed lines in Fig. partition the phase plane into 
regions where x and s have different sign combinations. In region I, the rate 
of biomass growth exceeds the rate of death (i > 0), and the rate of substrate 
inflow exceeds its rate of consumption and removal (s > 0); thus both biomass 
and substrate concentrations increase. However, eventually the biomass concen- 
tration becomes high enough that the rate of substrate consumption (combined 
with its outflow) exceeds its rate of inflow; the system then enters region II where 
the substrate concentration is decreasing (s < 0), but remains high enough that 
biomass continues to accumulate {x > 0). When the substrate concentration 
decreases to the point that the biomass growth rate is less than its net loss 
rate, the system enters region III: here the biomass concentration decreases 
{x < 0), and the substrate concentration continues to decrease (s < 0). When 
the biomass concentration has decreased sufficiently, however, the rate of sub- 
strate consumption becomes low enough that the substrate inflow rate exceeds 
its net loss rate, and the substrate concentration starts to increase (s > 0), while 
the biomass concentration continues to decrease (i < 0); this corresponds to 
region IV in Fig. Finally, the increase in substrate concentration causes the 
the biomass growth rate to increase, such that eventually the biomass growth 
rate exceeds its net loss rate, and the system again enters region I. This oscil- 
lation mechanism relies on the biomass relaxing more quickly or on a similar 
timescale to the substrate concentration (X^ > A^); rapid changes in biomass 
concentration lead to deficits or excesses in the substrate concentration, which 
are only slowly restored towards equilibrium by the flow of substrate into or out 
of the system. 

In contrast. Fig. |4j3 shows phase plane trajectories for a case where our 
analysis does not predict oscillations (/? = 10, x = 0). In this case, when /? -|- 
X ^ 1 the substrate concentration relaxes much more quickly than the biomass 
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Figure 4: Phase plane plots of the dynamics of microbial growth for a small substrate influx 
rate (/3 = 0.1) that gives rise to damped oscillations (A) and for a large substrate influx rate 
(/3 = 10) that gives exponential relaxation (B) to steady state (in both cases X = 0). In panel 
A (/3 = 0.1) oscillations arise due to a cyclical imbalance of growth vs death and net inflow 
vs consumption. In contrast, in panel B (/3 = 10) the substrate relaxes quickly to give a 
quasi-equilibrium of inflow vs consumption; biomass and substrate subsequently evolve slowly 
to the fixed point. 



concentration (A.; 3> A^;). Starting from given initial substrate and biomass 
concentrations, first the substrate concentration adjusts such that its inflow and 
net loss rates are equal (i.e. s « 0), for the given biomass concentration, then 
the biomass concentration slowly relaxes to the fixed point (making the quasi- 
equilibrium assumption that 6s' ~ 0, we find dx ~ —X^Sx), with accompanying 
changes in the substrate concentration. In other words, the trajectories rapidly 
approach the nullcline s — 0, then more slowly move along this nullcline to the 
fixed point. 



3. A single population: stochastic model 

We next consider the dynamics for the same model, but for small microbial 
populations. In this case, randomness in the birth and death/removal of indi- 
vidual microbes gives rise to stochastic fluctuations (or "demographic noise"), 
which cannot be neglected. McKane and Newman ^] have shown that, for a 
simple two-species predator-prey ecosystem, this demographic noise can produce 
sustained stochastic oscillations where the equivalent deterministic dynamical 
system shows only transient oscillations. Here, we show that the same effect 
happens in our single-population model. 

Our system can be described by the following chemical reaction scheme: 

0^5 (13) 

X + jS-^2X (14) 

X ^ (15) 

5^0 (16) 
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in which X and S denote microbes and substrate molecules respectively; we 
denote the numbers of microbes and substrate molecules as nx and ns re- 
spectively. Substrate molecules enter the system at a constant rate b' (Eq fTS)) . 
Microbes replicate upon consuming 7 molecules of substrate (EqlT4| , at a rate 
given by the Monod function fJ.'{ns)nx = T^rj^'fT-x ■ In addition, microbes 
are removed from the system at rate Dnx fEa llSI) and substrate molecules are 
removed at rate Rns (Eq \W\\ . The new parameters b' and K' arise because 
our units are now absolute numbers of microbes and substrate molecules rather 
than concentrations as in Eqs.(IT|) and ([2]); these parameters can easily be re- 
lated to b and K in our deterministic model 0- The parameters v, R, D and 
7, which do not depend on the units of substrate and biomass, are the same as 
in our deterministic model. We assume that reactions (|13p to are Poisson 
processes - i.e. they happen randomly in time with average rates determined by 
their rate constants. 

We use stochastic simulations to generate dynamical trajectories correspond- 
ing to the reaction scheme ([T3| - ([T6|) . In principle, this could be done using a 
kinetic Monte Carlo scheme such as the Gillespie algorithm [4l| , in which a single 
reaction happens in each timestep. However this would be extremely inefficient 
because the number of substrate molecules is typically very much larger than 
the number of microbes, so that reactions ([T3| and ([T6l) would happen much 
more often than reactions (fT4|) and (fT5|) . Instead, we map the system onto a set 
of differential equations for the concentrations of biomass and substrate, as in 
Section [21 but including a stochastic noise term which accounts for the effects 
of fluctuations (i.e. a set of Langevin equations for the biomass and substrate 
concentrations) . This is done using a Kramers-Moyal expansion of the Master 
Equation corresponding to Eqs (fT3| - (fT6| , to obtain a non-linear Fokker-Planck 
equation, which is an approximate description of the stochastic dynamics. We 
then write down the equivalent non-linear Langevin equation whose stochastic 



trajectories correspond to this Fokker-Planck equation [41|, |42|. Carrying out 
this procedure, as detailed in [Appendix B we find 



^=A(0) + Bi/2(0)^(t) (17) 

where cf) = {x, s)^ is the vector of concentrations and the noise vector ^(t) is 
Gaussian with zero mean and moment (^^^) — \6{t — t')/V, where I is the 
identity matrix. The vector A, given by 



xfi{s) — Dx 
—^xfi(s) + b — Rs 



(18) 



describes the deterministic time evolution of the system, as in Eqs.(lT]) and © 
(/x(s) = -j^^ being the Monod growth function). The matrix B, which ensures 



^For example b' (in moles per hour) is given by bV where b is measured in moles per litre 
per hour and the volume of the system V is in litres. 
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the correct coupling between biomass and substrate fluctuations, is given by 




We generate trajectories corresponding to this Langevin equation using an Euler 
integration scheme. 



A 



B 




Time (hrs) 

Figure 5: Stochastic time series for the growth of Escherichia coli on glucose, obtained by 
numerical simulation of the Langevin Equations 1171 1. for parameters v = lhr~^, D = 0.5hr~^, 
K = 7 = 1.8 X 10^" substrate molecules. Panel A shows the biomass concentration, 

normalized by its average, while panel B shows the normalized substrate concentration. Re- 
sults are shown for R = and b = 0.01, 0.1, 1 and 10 /iMhr^^, corresponding to x = Oi 
/3 = {0.01,0.1, 1, 10} (these parameter combinations are shown as green circles in the /9 — x 
phase diagram of Figure [!}. The insets show the fluctuations of biomass (A) and substrate 
(B) for /3 = 1 and /3 = 10, in more detail and on a shorter timescale. The colour codes in 
panels A and B are the same. 

Figures [5] & |6] show the resulting stochastic dynamical trajectories, for the 
same parameter set used for the deterministic trajectories of Figure [H repre- 
senting Escherichia coli growing on glucose, for a system volume V = 1ml. 
Trajectories are plotted for combinations of 6, R and D corresponding to the 
green circles in the phase diagram of Figure [TJ The average number of sub- 
strate molecules represented by each of these simulations is Inmol and the 
average number of microbes varies from sa 60,000 at 6 = 10/iM/hr to « 60 at 
h = 0.01/iM/hr (taking R = 0): the microbial population is indeed small enough 
(especially at low substrate inflow rates) that we would expect demographic fluc- 
tuations to play a significant role. Figures [5] & [6] show the microbial biomass 
and substrate concentrations (panels A and B respectively), normalised by their 
steady-state time average. For the two parameter combinations which lie outside 
the region of predicted oscillations in Figure [T] - a high rate of substrate influx 
(x = 0, /3 = 10, Figure[5]) or rapid substrate removal rate (x = 1, /3 = 0.1, Figure 
[6]) - the stochastic simulations show random fluctuations about the steady state 
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Figure 6: Stochastic time series for the growth of Escherichia coli on glucose, obtained by 
numerical simulation of the Langevin Equations l|17|l . for parameters v = lhr~^, D = 0.5hr~^, 
K = 7 = 1.8 X 10^'' substrate molecules. Panel A shows the biomass concentration, 

normalized by its average, while panel B shows the normalized substrate concentration. Re- 
sults are shown for /3 = 0.1 and x = {0, 0.25, 1} (green circles in Figure [TJ. These parameter 
combinations correspond to 6 = 0.1/iM/hr, R = (red line), b = 0.225/^M/hr, R = 0.125hr"l 
(brown line) and b = 0.6/iM/hr, _R = 0.5hr~^ (green line) respectively. The colour codes in 
panels A and B are the same. 



with no tendency to oscillate. However, for all other parameter combinations, 
which lie inside the region of predicted oscillations in Figure [1] the stochastic 
simulations show sustained oscillations. Comparing the deterministic trajecto- 
ries of Figure [5] with the stochastic trajectories of Figure [S] & IHl we see that, at 
least for these parameter combinations, sustained oscillations in the stochastic 
system occur for parameter sets where the deterministic system shows transient 



oscillations. Following McKane and Newman [27[, we reason that these sus- 
tained stochastic oscillations are generated by the continuous excitation of the 
oscillatory modes of the system by the intrinsic demographic fluctuations. 

To analyse in more detail the nature of the oscillations, we plot in Figures 
[7] and [5] (squares) the power spectrum of biomass and substrate concentration 
fluctuations, for stochastic simulations with the same parameters sets as in Fig- 
ures [5] and [51 Assuming a stationary stochastic procesf0, the Wicner-Khinchin 
theorem [42], relates the power spectral density I{uj) to the Fourier Transform 
of the autocorrelation function of the concentration fluctuations in the steady 
state - for example for the biomass concentration 

/>oo 

Ix{uj)= {Sx{t)Sx{t + T))e''^^dT (20) 



where we define the concentration fluctuation Sx{t) relative to the average (x): 
6x{t) = {x{t) — (x)). An equivalent formula holds for the power spectral density 



^A stationary stochastic process is defined to be one whose properties do not change with 
time and for which the autocorrelation function only depends on the difference between two 
time points r and not on their absolute times. 
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of the substrate concentration fluctuations. A peak in the power spectral 
density at a given frequency is a signature of sustained oscillations at that 
frequency in a noisy time series. Figures [7] and [8] clearly show that sustained 
oscillations indeed arise for those parameter sets for which the deterministic 
model produces transient oscillations (see Figures [U [2] and |4|) . 

We can also obtain analytical predictions for the power spectral density, 
using the van Kampen system size expansion [42] of the Master Equation corre- 
sponding to Eqs (fT3|) to p6)) . This procedure, the details of which are given in 
[Appendix C[ results in the following expressions for the power spectral densities 
of the biomass and substrate oscillations: 

r A W + xr + x' + 2c.yD' 

and 

^'^""^ " {(3 - ujyo^y + [p + x?^VD^ ■ ^^^^ 

where A = 13{V ^s^J\s=s■')^^ ■ Eqs (pij) and (1^^ . which are plotted in Figures 
[7] and [5] (solid lines) are in excellent agreement with the simulation results 0. 
Note that the amplitude of the PSD varies as , which indicates that the size 
of the fluctuations decrease as the volume of the system increases, as expected. 
These analytical results also allow us to determine which combinations of /3 and 
X give rise to a peak in the power spectrum, and hence to sustained stochastic 
oscillations. The red line in Fig. [T] shows the region of the /3-x parameter space 
in which sustained stochastic oscillations are expected from Eqs ([2T|) and (f22|) . 
Interestingly, this region lies inside the region where transient oscillations are 
predicted for the deterministic model. This suggests that, at least for this 
system, the presence of transient oscillations in the deterministic model is a 
necessary, but not sufficient, condition, for the presence of sustained oscillations 
in the stochastic model. 

As well as predicting the presence or absence of stochastic oscillations for a 
given parameter set, it is also important to know how pronounced these oscil- 
lations are expected to be. For example, in Figure [Sj oscillations are present 
both for /3 = 0.1 (red line) and (3 = 0.01 (blue line), but they are much more 
pronounced in the latter case. Similarly, for (3 = 0.1 and /3 = 0.01 we see a 
corresponding increase in sharpness of the peaks in the power spectra. We can 
understand this in terms of changes in the quality factor of the oscillations ex- 
amined in Section [2l As (3 is decreased (for x — ^) the ratio of the timescales of 



^In Figure [T] the deviation at high frequency between the analytical and simulation results 
for the substrate power spectrum for small /3 (corresponding to low substrate inflow rate) may 
be attributed to the fact that here the relative fluctuations in biomass and substrate are both 
roughly of order ^ 1, so that one would expect the Kramers-Moyal expansion used to obtain 
the simulation algorithm to give different results from the van Kampen expansion used to 
obtain the expression for the power spectrum. 
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biomass and substrate relaxation decreases, leading to an increase in the qual- 
ity factor as predicted by Ea.(IT2|). More generally, as Fig|3]and Eq. lfTTt shows, 
the condition /3 + x <C must be satisfied to show significant oscillations; 
for example, in FiglU we see that for /3 = 0.1 the peak in the power gradually 
disappears as x is increased due to a decrease in the relaxation time of the 
substrate relative to biomass. 




lo" 10^ 10^ lo' 10° lo' 

Frequency (hrs~^) 



Figure 7: Simulation (points) and analytical (solid lines) results for the power spectrum of 
biomass (A) and substrate (B) concentration fluctuations, for substrate influxes h of 0.01, 0.1, 
1 and 10 /iMhr~^, corresponding to /3 = {0.01, 0.1, 1, 10}, with x = in all cases. Clear peaks 
in the power spectrum for fi = {0.01,0.1,1} indicate the presence of sustained stochastic 
oscillations. Note that the colour codes in panels A and B are the same. 



4. A two-species nutrient-cycling ecosystem 

Our results thus far demonstrate that demographic noise can result in sus- 
tained oscillations for a single microbial population whose rates of biomass and 
substrate removal are not balanced. In the natural environment, however, the 
waste product of one microbial population may form the substrate for another, 
leading to microbial ecosystems with complex webs of crossfeeding interactions. 
Would we expect to see oscillations in such complex ecosystems? 

Figure [H] shows a schematic illustration of a simple nutrient-cycling microbial 
ecosystem. Here, population 1 consumes substrate 1 and produces substrate 2; 
this is in turn consumed by a second population, which releases as its waste 
product substrate 1. We assume that only substrate 1 is supplied by the external 
environment. This model might represent for example the cycling of carbon 
between methane and carbon dioxide by methanogens and methanotrophs, or 
the cycling of sulphur between sulphide and sulphate by sulphur oxidising and 
sulphur reducing bacteria; here, the environment is assumed to supply one of 
the forms of carbon/sulphur at a constant rate, and other necessary inputs 
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Figure 8: Simulation (points) and analytical (solid lines) results for the power spectrum of 
biomass (A) and substrate (B) concentration fluctuations, for /3 = 0.1 and x = {0,0.25,1}, 
corresponding to the (6, J?) combinations (0.1/iM/hr, 0), (0.225/iM/hr, 0.125hr~l) and 
(0.6^M/hr, 0.5hr~^), respectively. Clear peaks in the power spectrum for x = {0,0.25} 
indicate the presence of sustained stochastic oscillations. Note that the colour codes in panels 
A and B arc the same. 



Inflow of 
Substrate 1 




Figure 9: Schematic illustration of a two-species nutrient-cycling ecosystem. 



(e.g. oxygen for the oxidation reaction and hydrogen for the reduction reaction) 
are assumed to be available in plentiful supply. The cycle of Figure can be 
represented by the following set of deterministic equations: 
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dt 

dX2 

'dt 
dsi 
"dT 
ds2 
'dt 

Here, Xi and X2 are the biomass concentrations of populations 1 and 2, while 
Si and S2 are the concentrations of the two substrates. The factor / arises in 
the equations for dsi/dt and ds2/dt because we suppose that for every gamma 
substrate molecules which they consume, microbes use a fraction (1 — /) for 
growth and excrete the remaining / in the form of waste product, which ensures 
conservation of matter. We assume for simplicity that both populations have 
the same growth/removal parameters v, K, 7 and D, and that the substrate 
removal rate R is zero. The inflow rate of substrate 1 is given by b. Eqs. d^^ 
have a single non-trivial fixed point with steady state biomass and substrate 
concentrations given by 



b 




P) 


bf 






P) 


K 




{v/D~ 


1) 



(24) 



Fig. [To] (panel A) shows deterministic trajectories for the biomass concentra- 
tions of microbial populations 1 and 2 predicted by Eqs. (P^ . using the parame- 
ter set for Escherichia coli growing on glucose defined in Section[21 with / = 0.8. 
Results are shown for two different substrate inflow rates. As we observed for 
the single population in Figure [U transient oscillations occur for the lower sub- 
strate inflow rate b = 0.01/iM/hr; these disappear on increasing the substrate 
inflow rate to 6 = 10/zM/hr. We hypothesise that the underlying mechanism 
for these oscillations is the same as that for the single population model: i.e. 
transient imbalances in the biomass and substrate concentrations which arise 
because the biomass concentrations change on a comparable or faster timescale 
than the substrate concentrations. 

Using the same procedures as in Section [3] and [Appendix B[ we can also 
carry out stochastic simulations for this model. The resulting trajectories are 
shown in Fig. [TOBfcC. for a system volume oi V = 1ml and the same param- 
eters as for the deterministic case. Sustained stochastic oscillations are indeed 
observed for the substrate inflow rate b = O.Ol/xM/hr (Panel B) which arise 
from the transient oscillations in the deterministic model (Fig. ITOR - dashed 



K + si 
VX2S2 

K+S2 



-DX2 



VXlSl VX2S2 

'^K + si ■'"^K + S2 

VX2S2 . VXiSi 
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Figure 10: Stochastic population oscillations for the mutualistic cycle model (Fig. O with 
V = lhr-\ D = 0.5hr-\ K = l^lM, 7 = 1.8 x 10^'', / = 0.8 and volume V = 1ml. Blue lines 
correspond to xi and red to X2. Panel A shows deterministic times series for inflow rate of 
substrate 1 of 6 = O.OlAtM/hr (/3 = 0.01 - dashed Unes) and of 6 = lO^tM/hr (/3 = 10 - solid line) 
calculated by numerical integration of Ean l23l Panels B and C show the equivalent results, 
but from stochastic simulations, for b = 0.01/jM/hr and b = lO/xM/hr, respectively. We see 
that at the lower inflow rate (b = 0.01/iM/hr), damped oscillations (A) in the deterministic 
case become sustained stochastic oscillations in the stochastic model (C). However, for the 
higher inflow rate (b = 10/iM/hr) neither the deterministic (A) nor the stochastic (C) model 
shows oscillations. These results mirror the observation of stochastically forced oscillations at 
lower inflow rates in the single population model of Section \3\ Note that the colour codes in 
B and C are the same as in A. 



lines), while for the higher inflow rate b = 10/iM/hr, neither the stochastic 
(Fig. \WC) nor the deterministic system shows oscillations (Fig. [TOR - solid 
lines). Fig. [Til shows biomass concentration power spectra for the stochastic 
simulations: as expected, for b = 0.01/iM/hr a clear peak is present, indicat- 
ing sustained stochastic oscillations, while the power spectra show no peak for 
b = 10/L(M/hr. 

The nutrient-cycling model investigated here represents, of course, only one 
example of many possible microbial ecosystem topologies. Detailed analysis 
of the effects of ecosystem topology on the occurrence and character of noise- 
induced oscillations would be an interesting subject for further work; neverthe- 
less this example suffices to show that the phenomenon is likely to be widespread 
in microbial ecosystems with small population sizes and unbalanced biomass and 
substrate removal rates. 
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Figure 11: Power spectra of the stochastic time series shown in Fig. llOC fcD for the mutuahstic 
cycle model. Blue squares correspond to xi and red to X2', filled symbols correspond to 
substrate inflow rate b = O.Ol/^M/hr, while open symbols correspond to 6 = 10/xM/hr. The 
presence of a peak in the biomass power spectrum for b = 0.01/iM/hr indicates sustained 
oscillations; no such peak is present for the higher influx rate b = 10/iM/hr. 



5. Discussion 

In this paper, we have demonstrated that microbial populations can show 
oscillatory dynamics, under conditions where the biomass is removed from the 
system faster than the substrate, and the supply rate of substrate is not too 
high. For large populations these oscillations are transient, and the population 
eventually reaches a non-oscillatory steady state. For small populations, how- 
ever, the oscillations can be sustained indefinitely. The origin of these sustained 
stochastic oscillations is the mechanism discovered by McKane and Newman 
[23], in which underlying oscillatory modes of the deterministic system are ex- 
cited by intrinsic demographic noise. While this mechanism has been shown to 
produce stochastic oscillations in a number of different systems [13-01, to our 
knowledge this is the first demonstration of its relevance in a simple model of a 
single population growing on a substrate. For the single population model, we 
have used both simulations and analytical arguments to investigate in detail the 
conditions for the presence of these oscillations and their characteristics. We 
have also shown using simulations that the same effect can be expected to arise 
in more complex microbial ecosystems with cross-feeding interactions. 

What are the implications of these results for microbial populations in the 
natural environmnent? Our work suggests that oscillatory dynamics may be 
widespread in circumstances where microbial populations are small, with high 
death rates. Small microbial populations are likely to be found in microenviron- 
ments such as the interstices between soil grains, on particulate matter in the 
ocean or, for pathogenic bacteria, the inside of a host cell. In such a closed envi- 
ronment, it is indeed likely that the removal rate of substrate may be low, while 
the death rate of microbes due to phage predation or attack by host defence 
mechanisms, may be high. If such oscillations are indeed widespread, this would 
have important implications for the establishment and maintenance of microbial 
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communities: for example, oscillations would be likely to have a strong effect 
on the distribution of times to extinction. There may also be an interesting 
interplay between microbial ecosystem topology and the characteristics of the 
stochastic oscillations. 

The next stage in this work is clearly to test our predictions experimentally. 
Traditionally, theories about microbial population dynamics are tested using 
chemostats, but in the absence of mortality factors such as phage the conditions 
required for oscillations are not satisfied in a conventional chemostat, and typi- 
cal population sizes are anyway almost certainly too large to observe sustained 
noise-induced oscillations. Microchemostats however, in which small microbial 
populations are maintained under steady-state conditions in microfluidic de- 
vices (43| . may provide conditions under which noise-induced oscillations could 
be observed. Moreover, microscopic techniques for observation of the growth 
of microbes in complex confined geometries 44] or as biofilm communities [Zsj ] , 
provide exciting possibilities for testing the likely significance of stochastic os- 
cillations for microbial communities in the natural environment. 



Acknowledgments 

The authors thank Richard Blythe and Mike Gates for valuable discussions 
and reading of the manuscript, and Billy K. Huang for his contribution to 
the early stages of this project. This work was funded by the Leverhulme 
Trust under grant number F/00158/BX, and by EPSRC under grant number 
EP/E030173. RJA was funded by the Royal Society of Edinburgh and by a 
Royal Society University Research Fellowship. AF thanks Prof. Kenneth Mur- 
ray and the Darwin Trust for fellowship and additional support. 



Appendix A. Quality factor of the biomass oscillations 

The quality factor of a resonant system measures the fidelity of its oscilla- 
tions; for a mechanical or electrical oscillator, the Q factor can be interpreted 
as the ratio of the energy stored to energy dissipated over an oscillation cycle 
(iij . When the Q-factor is large an approximate definition is given in terms 
of its power spectrum /(w) as the ratio of the peak oscillation frequency rip to 
the range of frequencies that significant oscillations occur Aw: Q « Jlo/Aw. A 
simple damped harmonic oscillator such as a mass on a spring in a viscous fluid, 
which oscillates at frequency with no friction, has the equation of motion: 

d^x C dx n ... 
— + ^—+ujIx^Q A.l 
at'^ m at 

where m is the mass, x is its position and Q is the friction coefficient. For the 
damped oscillator, the quality factor Q is proportional to the ratio of energy 
stored in the spring, to the energy dissipated to frictional loss, per cycle: 

Q « ^ (A.2) 
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For the problem under consideration in this paper, our hnear approximation of 
the deterministic dynamics, Eqs ([5]) and ([6]), aUows us to write 



— w Ss A.3 

at 7 

— ^ « -D-fSx ~ D{I3 + x)5s 
at 

which can be combined into a single second-order differential equation for the 
biomass concentration close to the steady state: 

^ + + X) ^ + D^pSx ^ (A.4) 

Comparing Eqs (lA.ll) and (jA.4[) we see that our system maps onto the damped 
simple harmonic oscillator, with D(Ji + x) and playing the roles of C^jm and 
respectively. We can therefore use this analogy to predict that the quality 
factor will be given by 

As shown in the main text, the quality factor can be expressed as Q = 



^JX^jXa-, such that the requirement for high Q oscillations is >• A.;. Intu- 
itively, (for small fluctuations around the fixed point) the relaxation rate of the 
substrate A^ = /3-|-x is analogous to the rate that a mechanical system dissipates 
energy C/™j and so we can understand the presence or absence of oscillations 
of biomass and substrate as analogous to the underdamped and overdamped 
limits of a SHO. 



Appendix B. Langevin approximation of the stochastic dynamics 

In this Appendix, we briefly describe the derivation and implementation of 
the Langevin approximation of the stochastic model, via a Kramers-Moyal ex- 
pansion, which we use in our simulations. We follow the procedure described 
by Van Kampen [i^l • We begin by writing the Master Equation for the proba- 
bility P{nx ,ns,t) of observing the system with nx microbes and ns substrate 
molecules at time t. Denoting {nx,ns) as the vector n, the Master Equation is 
given by 

= ^ a',{n + ru)P{n -f r^, t) - a'^{n)P{n, t) (B.l) 

fc=i 

In Eq. (jB.l[) . the stoichiometry vector r^j = {rkXifks) denotes the change in n 
when reaction k flres: ri = (0, 1), r2 = (1, —7), — (—1, 0) and r4 — (0, —1). 
The propensity function a'^, denotes the probability of occurrence of reaction k 
per unit time: a'l = 6', = ^'{ns)nx^ a'^ = Dnx and 04 = Rns- The flrst 
term in Eq. ()B.l[) represents the flux of probability from other states n + 
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into state n, while the second term denotes flux of probabihty due to reactions 
which move the system out of state n. 

We now Taylor expand a'i^{n + rk)P{n + rfc, t) about n, to second order in 
r/c, and substitute the result into Eq. (jB.l[) : this leads to the nonlinear Fokker- 
Planck equation 



dP{n,t) _ ^ 
fe=i 



dt 



rkB 



dqk 
dnx 



k=l 



kX 



d'^nx 



2rkxrks 



dns 
d^qk 



(B.2) 



dnxdns 



' kS 



d^qk 



where qk = a'f.{n)P(n, t). This approximation will be good as long as the num- 
ber of microbes and substrate molecules is large (see comments in |Appendix C[ ) . 
This in turn can be represented by an equivalent nonlinear stochastic differential 
equation (Langevin Equation), for which probability distribution of trajectories 
generated will follow Eq. (|B.2|) . This Langevin Equation is given by [42| : 



dn 

d^ 



A'(n) + B'i/2(n)|'(t) 



(B.3) 



where the vector A' is given by 



A' 



E 



the matrix B' is given by 



B' 



^a'^rkrl. 

k=l 



(B.4) 



(B.5) 



and ^' (t) is a vector of Gaussian (white) noise components, with zero mean and 
second moment (^'^'"^) = \S{t—t'), where I is the identity matrix. Note that as B 
is a symmetric matrix, V is unitary (i.e. = V~^), which implies B^/^ — B^^^^ . 
This means multiplying B^/^ by an arbitrary unitary matrix leaves the variance 
matrix B, and hence, the probability distribution of paths is unchanged, as the 
Fokker-Planck equation only depends on B. We can transform this Langevin 
equation in terms of the stochastic dynamics of the concentrations, x = nx/V 
and s = ns/V, by dividing through by the volume of the system V. The result 
is that the propensities Ofc = a'f./V, so that oi = 6, 02 = n{s)x, 03 = Dx and 
04 = Rs, leading to the Langevin equation in the main text fEgnlTT)). with A 
and B given by Eqns [TH] and [12] and where the noise term ^ = C/W. 

In our simulations, we compute at each timestep the propensities for 
each reaction k, and hence obtain the vector A and matrix B. We then compute 
the square root of B numerically, by calculating the eigenvalue decomposition 
of B = VDV''', to give B^/^ = VD^/^N/t, where D is the diagonal matrix of 
eigenvalues and V is a matrix whose columns are the eigenvectors of B and f is 
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the transpose and conjugate operation (Hermitian conjugate). The components 
of the vector ^ are obtained using a random number generator. This allows 
us to update the microbe and substrate populations using an Euler integration 
scheme with timestep of 0.05hrs. 



Appendix C. Analytical calculation of the power spectrum 

To obtain analytical expressions for the power spectrum of biomass and sub- 
strate concentration fluctuations, we need to find a way of expressing dynamics 
of the fluctuations around the fixed point. To do this we apply van Kampen's 
large system expansion method [i^l to the Master Equation (IB.ip . This method 
relies on the ansatz that the solution of Eq. (|B.H) will be of the form 



n{t) = V<jyo{t) + WiPit) (C.l) 

The vector function 4'Q{t) is the solution of the deterministic equations ([1} and 
or equivalently 

d0o/di = A, (C.2) 

while the vector function xp{t) represents the fluctuations around the determin- 
istic solution, caused by noise, which we expect to be of size ^/V. The Van Kam- 
pen method then formally expands the Master eqn in orders of V^^^^ about the 
deterministic solution to give a linear noise Fokker Planck or Langevin equation. 
We can do this by starting with the Langevin equation from the Kramers-Moyal 
expansion fEan lB.3l) and substituting in Egn lC.fl on both sides of the equation 
to give: 

+ B"/=(V0„ + •/Vi,)('{t) 

Each of the A' and B' are composed of terms a'^r^., a Taylor expansion of which 
gives 



(C.4) 



where the contribution to the Jacobian of each reaction is given by = 
rfe(Va5j)-^; these are related to the Jacobian by J = J2k^k- Using Eqn lC.4l 
in Eqn IC.Sl and keeping only terms of order V and VV and assuming that V 
is large so we only retain terms linear in the noise (known as the Linear Noise 
Approximation), we have 

^BV2(0^)|'(t), 
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where we have used the fact that a'f^{Vcf)Q) = Vak{4>Q). Further, using EQn lC.2l 
means we can cancel the terms of order V that relate to the deterministic or 
'macroscopic' evolution of the system to leave a Langevin Equation for the 
stochastic part: 

^^^^P + B'/\4>,)^'{t). (C.6) 

However, we are interested in the fluctuations in concentration about the steady 
state or fixed point, so we can form the variable z — ip/^/V — {Sx,Ss)^ and 
consider the solution for (f>Q(t — > oo) = (p* to give an effective linear-noise 
Langevin equation for the fluctuations around the steady state: 

d z 

— = J*2 + B*i/2^(t) (C.7) 

where J* is the Jacobian matrix of the dynamics around the fixed point, given 
in Eqniniand the variance matrix B* — B(</)q) of the effective diffusion process 
is Egn lB . 51 evaluated at the fixed point: 



B* =^Afcrfcr 
fe=i 

= (6 - Rs*) 



(C.8) 

2/7 -1 



b-Rs' ^ I 



The Langevin noise term is Gaussian with zero mean and moment (^C^) = 
1(5 (t — t')/V. The van Kampen and Kramers-Moyal approaches produce linear 
and non-linear Fokker-Plank or stochastic differential equations. However, as 
we have shown they can be obtained from each other by a change of variable 
and are equivalent to within fluctuations of order VV, which in any case is the 
level of approximation of each approach. 

The power spectrum of fluctuations around the steady state solution can be 
calculated from the Fourier Transform (FT) of Eqn lC.7l as 

(Z(cj)Z(w)t) = ^Kiu)B*K\Lu) (C.9) 

where Z(a;) = FT {z{t)} and K{uj) — (iwl — J)~^. In particular, if we assume that 
7^1 then we flnd the power spectrum of biomass and substrate fluctuations 
are as given by Eqns. [2T]&[22]in the main text. 
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